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We show, using an ah initio approach based on Quantum Monte Carlo technique, that the pseu- 
dogap regime emerges in ultracold Fermi gases close to the unitary point. We locate the onset of 
this regime at a value of the interaction strength corresponding to ~ —0.05 (a - scattering 

length) . We determine the evolution of the gap as a function of temperature and interaction strength 
in the Fermi gas around the unitary limit and show that our results exhibit a remarkable agreement 
with the recent wave- vector resolved radio frequency spectroscopy data. Our results indicate that a 
finite temperature structure of the Fermi gas around unitarity is more complicated and involves the 
presence of the phase with preformed Cooper pairs, which however do not contribute to the long 
range order. 
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The unitary Fermi gas is a dilute but an exceptionally strongly correlated ensemble of particles exhibiting universal 
properties and also a superfluid with a very large critical temperature, making it relevant to a wide range of systems 
including the quark-gluon plasma, neutron stars, nuclei [1], and to a certain extent to high Tc superconductors. A 
number of fundamental properties of the Fermi gas close to the unitary limit remain unknown and still pose an 
unprecedented challenge for the theory [2] . Although a convincing experimental proof of superfluidity was provided 6 
years ago [3] , the measurements of the pairing gap at finite temperatures are still in their infancy and only recently the 
technique of photoemission spectroscopy allowed to probe the low-energy excitations of the gas at finite temperatures 
Q . The pairing gap measurements are believed to answer one of the most intriguing questions related to the strongly 
interacting Fermi gas: Does the simple picture of the phase transition from superfluid to normal state still holds 
around the unitary regime? Here we show, using an ab initio approach based on Quantum Monte Carlo (QMC) 
technique, that the pseudogap regime emerges in ultracold Fermi gases close to the unitary point. In the pseudogap 
regime fermions still have a strong tendency to be correlated with each other and behave to some extent as a 
system of bosons. Consequently above the critical temperature the system exists in an exotic state which is neither 
entirely normal nor superfluid and thus clearly eludes the understanding within the framework of classical BCS 
theory of superconductivity. This situation can be better understood from the viewpoint of the deep Bose-Einstein 
Condensation (BEC) limit, where the pairs of fermions form bound bosonic dimers. This implies the existence of 
two distinct temperature scales: the first one, T c , related to the Bose-Einstein condensation of bosonic molecules and 
the second one, T* > T c , which is a temperature needed to break up a dimer into fermions. Somewhere around the 
unitary regime these two temperatures become practically indistinguishable and eventually merge into the critical 
temperature for superfluid to normal phase transition [5j. The region between T c and T* is usually referred to as 
a pseudogap phase, where the system, although not being a superfluid, still exhibits the gap in the spectrum of 
quasiparticle excitations. Actually the term "pseudogap phase" may be somewhat misleading as, first, it is a strong 
depletion of the density of states in the spectrum and second, it is strictly speaking a particular regime, not a phase, 
as there is no phase transition between pseudogap regime and the normal Fermi gas. 

The situation in the ultracold atomic gases at the unitary regime is somewhat similar to that encountered in 
high temperature superconductors (HTSC) [6[, although the pseudogap regime in the cuprates is not exclusively 
related to Cooper-pair formation [7]. To appreciate these similarities one has to consider the ratio between the 
critical temperature and the pairing gap at zero temperature. According to the celebrated BCS theory this ratio is 
approximately equal to 1.76 (for the temperature expressed in energy units, ks = 1). Any deviation from this value 
is an indication of the non-BCS superfluidity in the system. The origin of the deviation is usually related to the 
interaction strength and manifests itself in a rather large ratio of the pairing gap A to Fermi energy sf- It stands 
in contradiction to the implicit assumption of the BCS theory that the interaction is weak which imply A/sp <C 1. 
Clearly both ultracold gases and HTSC violate this assumption, see Fig. [U Morevover the cold fermionic atoms 
constitute the system with the largest ratio A/sp and known HTSC only approach this limit |8j- 

Theoretically the single particle gap in the fermionic spectrum can be determined from the spectral weight function 
A(p, cj), which carry information about the single particle excitation spectrum. In our approach we have used the 
Path Integral Monte Carlo (PIMC) technique on the lattice, which provides numerically essentially exact results with 
quantifiable uncertainties. The details of the approach have been described elsewhere [9l4l2j. Calculations presented 
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FIG. 1: (Color online) Ratio of the pairing gap and the critical temperature for various superconductors (data extracted from 
Ref. [3(3]) . The clearly visible deviation from the typical ratio predicted by the BCS theory is a hallmark of HTSC. The same 
ratio for atomic Fermi gases close to the unitary limit are denoted by blue circles (data from Ref. |9|]). 



here have been performed on a larger lattice than reported in Ref. [11], namely 10 3 , with particle number varying 
between 90 and 110. The number of Monte Carlo samples was chosen in such a way to decrease the statistical errors 
below 1%. The systematic errors, some due to finite lattice effects, others due to finite range effects, are estimated 
at no more than 10%, depending on physical quantity. In order to determine the size of the gap, we have used the 
one-particle temperature Green's (Matsubara) function which can be easily calculated within the PIMC formalism 
and which is related to the spectral weight function [13]. The PIMC method provides us with a discrete set of values 
of the Green's function. In the present calculations we have determined the Green's function at 51 imaginary time 
points, equally spaced within the interval [0, j3 = 1/T]. We have found that due to the smooth behavior of the Green's 
function the accuracy is not affected when increasing the number of points. We have used two methods to extract 
the spectral weight function from the Matsubara function: the maximum entropy method (MEM) and the Singular 
Value Decomposition method (SVD) [Tll.[l^. 
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FIG. 2: (Color online) Spectral weight function A(p,oj) at l/kpa = 0.2 for three temperatures: T = 0.13sf < T c (left panel), 
T = 0.19sf ~ T c (middle panel) and T = 0.26s f > T* (right panel). White circles denote the position of the maxima of 
two branches of spectral weight function and the white line represents the quasiparticle energy E(p) fit to the maxima of the 
spectral weight function. In the insets we show the section of the spectral weight function at the Fermi level. 

As an example in the Fig. [2] the spectral weight function is shown on the BEC side of the unitary regime for 
(kpo) -1 = 0.2. The three temperature values correspond to three regimes realized in this system. At T w 0.13£f < 
T c = 0.19(l)£i? the system is superfluid [9], at T & 0.19s f the long range order is lost, but the gap in the single 
particle spectrum is still present. Finally, the lower panel shows the spectral weight function just above the crossover 
to the normal Fermi liquid. The negative energy branch is nonzero essentially up to p ~ pp. For the momenta above 
the Fermi surface the negative energy branch has a tail which in the case of the unitary Fermi gas is related to the 
behavior of the occupation probability n(p) at large p which decay as 1/p 4 [15]. The positive energy branch starts 
essentially around p ~ pp and the presence of the gap can be easily detected. The upper panel of Fig. [2] describes 
the superfluid system with no states at the Fermi surface (see inset). In the middle panel (T « T c ) the density of 
states at the Fermi level is significantly lowered. At large temperatures the two branches merge into a single branch 
extending from negative to positive values of lj — /i, characteristic for a normal state. Similarly like in [Til, llfij we 
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TABLE I: Pairing gaps A, self-energy U and effective masses m* extracted from the spectral weight function for various 
temperatures T and scattering lengths (see left column). In the left column the extracted values of crossover temperature T* 
have been shown together with the values of the critical temperatures (upper estimates) taken from Ref. [9|. All the results 
except of those for (/ci?a) _1 = —0.2 have been obtained for the lattice 10 3 . The results on the BCS side at (/ci?a) _1 = —0.2 
have been obtained for the lattice 8 3 as described in [Tl|. The values of A/sf < 0.1 should be treated as lying within the 
interval [0,0.1] due the finite resolution, see Ref. 0. 



have fitted the formula for the quasiparticle energy E(p) = yj (p 2 /2m* + U - nf + A 2 to the maxima of the spectral 
weight function. The extracted values of gap A, self-energy U and effective mass m* have been presented in the 
table HI One may notice that the effective mass is practically constant and does not deviate from the bare mass by 
more than 10%. The self-energy U decreases with (fc^a) -1 . When compared to the experimentally extracted value 
at unit arity (see Ref. |17|) it agrees within a 10% accuracy. The values of T c around unitarity have been determined 
elsewhere [9|, LL8j. The value of the gap at T c is usually called the pseudogap, i.e. measure of the strength of the short 
range correlations at the critical temperature. These results can be subsequently used to estimate the value of T*, 
where the gap tends to zero and also the value of the gap at the critical temperature T c for the superfluid to normal 
phase transition. We show the size of the gap A/sp plotted as a function both temperature and scattering length in 
Fig. [U The two temperatures T* and T c become indistinguishable, due to uncertainties generated by Monte Carlo 
approach, between —0.1 < (kpo) -1 < 0. Hence already at the unitary limit the system exhibits the pseudogap phase, 
albeit in a relatively small temperature interval. The temperature interval within which the pseudogap phase exists 
increases as one moves to the BEC side. 

The present work is the first attempt to determine the evolution of the pseudogap as a function of temperature 
and interaction strength from fully ab initio calculations. Previous results obtained using this ab initio method were 
used as a benchmark by other theoretical approaches and were as well validated by subsequent experimental results 
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FIG. 3: (Color online) The gap A/sf extracted from the spectral weight function as a function of temperature and scattering 
length. The dashed lines denote two temperatures: critical temperature T c and the crossover temperature T*. Uncertainties 
(both systematic and statistic errors, estimated to be no more than 10%) of these temperatures, are denoted by shaded area. 




FIG. 4: (Color online) Experimental (dots) and theoretical (lines) energy distribution curves (EDC) for the trapped atomic 
gas in the unitary regime at the critical temperature (defined in the center of the trap) for various momenta p. 



[l9|,[20|. On one hand, the most recent calculations of Ref. [2lj |. based on the pairing- fluctuation theory, indicated the 
presence of the pseudogap without however specifying its size and evolution as a function of the scattering length. On 
the other hand the recent measurements of the magnetic susceptibility of the trapped atomic gases were explained 
within the Fermi liquid theory without invoking the pseudogap model 2(J 22], a description challenged in Ref. 21]. 
The Fermi liquid parameters in Ref. [22| were extracted in a QMC calculation at zero temperature, in which the 
formation of pairs was artificially suppressed, and the spectral functions were approximated with (5-functions. Here 



we present the comparison of the experimental results based on the wave- vector resolved spectroscopy [21|, |23|, |2 



Since the experimental data are for a trap which is inhomogeneous and likely composed of spatially separated and 
coexisting phases it is difficult to relate the particular behavior of the extracted energy distribution curves to the 
appearance of the pseudogap phase. However the agreement of the experimental and theoretical results supporting 
the existence of the pseudogap provides an implicit confirmation of the quality of the theoretical calculations. We 
use the Local Density Approximation [25] to relate the PIMC input to trap data, see Fig. 2] (without using any 
fitting parameters, see 14] for details). Over the years there was a significant body of work addressing the physics of 
preformed pairs, pseudogap and dissociation temperature T*, see Refs. fBL l26l-[29| and earlier references, based on a 
variety of theoretical models and approximations. The theoretical predictions for some or all of these entities either 
varied widely or were at most of qualitative nature and the values obtained in these references for the pseudogap and 
T* are significantly different from ours. The lack of quantifiable error bars makes essentially impossible to re-conciliate 
these earlier predictions and to judge their accuracy. 

In summary we have determined in ab initio calculations the spectral weight function for various temperatures 
around the unitary regime. We have confirmed the existence of the pseudogap phase (we have reported the evidence 
of this phase at unitarity in Ref. [11]) and located the onset of this phase at a value of the interaction strength 
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corresponding to (fc^a) -1 « —0.05. To date the present results are the first ab initio calculations which were able 
to determine the size of the pseudogap and its behavior as a function of interaction strength and temperature. Our 
results indicate that a finite temperature structure of the Fermi gas around unitarity is more complicated and involves 
the presence of the phase with preformed Cooper pairs, which however do not contribute to the long range order. 
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SPECTRAL WEIGHT FUNCTION AND SINGLE PARTICLE PARAMETERS 



To extract the spectral weight function we have used the one-body temperature Green's (Matsubara) function which 
can be calculated within the PIMC formalism |l3j: 

GiP, r) = -|Tr{exp[-(/3 - r)(H - fiN)]^(p) exp[-r(H - fiN)]ft(p)}. (1) 

In the above expression j3 = 1/T is the inverse temperature and < r < /3. For the spin-symmetric system with 
the spin-independent interaction Q(p, r) is diagonal in the spin space and contributions coming from the spin- up and 
spin-down fermions are equal. Therefore we have suppressed the spin indices in all formulas. The trace Tr is performed 
over the Fock space, and Z = Tr{exp[— f3(H — pN)]}. The spectral weight function A(p,u) = + 
is defined as: 

A( + \p,lu) = — 5^exp[-^(£7 n -/iJV n )]<n|V'(p)|m>HV' t (p)|n)5(a; + E n - E m + fi), 

n,m 

A ( -\p,oj) = ^2exp[-^(£7 n -A^„)]<n|^t(p)| m )( m |^(p)|„)J( w + f; m -£7 n + ^) ) (2) 

n,m 

where \n) represents the state with energy E n and particle number N n . It is related to the temperature Green's 
function: 

G(p,r) = r du,A(p, U ) ** P( ~r\v (3) 
2tt 1 + exp(-w/3) 

and by definition A(p, u>) fulfills the following constraints: 

. , . f°° doj A , , , [°° dui A , 1 

A(p,u,)>0, j_^-A M = l, y_^_^ (P)W ) I _ w = n( p) I (4) 

where n(p) is the occupation number of a state with momentum p. The PIMC method provides us with a discrete set 
of values of the Green's function. In the present calculations we have determined the Green's function at 51 imaginary 
time steps, equally spaced within the interval [0,/3]. We have found that due to the smooth behavior of the Green's 
function the accuracy is not affected when increasing the number of points. The discretization of the left hand side of 
Eq. ((3]) leads to the class of numerically ill-posed inverse problems. We have used two methods, designed particularly 
to deal with such problems: the maximum entropy method (MEM) and the Singular Value Decomposition method 
(SVD). 

In the maximum entropy method [31(, which is based on Bayes' theorem, we treat the values of the Green's function 
calculated within PIMC approach as normally distributed random numbers G(p,Ti) (i — 0, 1,2, ...,A/* T = 50) around 
the true values 7$). The Bayesian strategy of searching for the most probable solution leads to the minimization 
of the quantity \x 2 ~ &S(M), where a > 0, 



x 2 



Y,[g(p,Ti)-g(p,Tii\ /a 2 , (5) 



i=l 

and S(M) is the relative information entropy with respect to the assumed model M: 

f A(p,u) k ) 

k 



S(M) = -^Auj A(p,uj k ) - M(u k ) - A(p,uj k )\n (■ 
h L V 



M(u k ) 



(6) 



The minimization is performed with respect to the values of A(p,Uk). The entropy term prevents excessive inclusion 
of unjustified structure into the shape of the spectral weight function. The constraints @ are enforced by means of 
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Lagrange multipliers. In our previous calculations obtained for the smaller lattice of size 8 3 , which were reported in 
Ref. [ll|, the calculated Green's function was not as smooth as the one we obtained for 10 3 lattice. Therefore in the 
present calculations we could decrease the value of a by an order of magnitude as compared to the previously used 
value: a 2 a « 0.3. It has substantially improved the resolution of the present calculations. We have used the energy 
window —5 < (lj — /i) / 'sf < 5 and calculated A(p,u)k) for 120 points equally spaced within this interval, although the 
size of the window and the number of points could have been chosen smaller without any effect on the results. In our 
MEM approach we have improved the method of finding the spectral weight function by constructing a sequence of 
minimizations with a gradually refined model. Namely, at the end of each minimization process the result has been 
used to defined a new model in the form of two gaussians which maximize the overlap: ^ k \/ A(p, Uk)M(uJk)Auj/27r. 
Note that since the spectral weight function and the model are both nonnegative and normalized, the above quantity 
can take values from the interval [0, 1] only. The parameters of the model A4 include the heights, widths and the 
relative shift of the gaussians. These parameters were adjusted to maximize the overlap and thus define a new model. 
Then the new minimization process was started. The procedure has been continued until the value of x 2 did not 
improve. The procedure has generated spectral weight functions, where both the position of the maxima and the 
widths were surprisingly stable, no matter what particular model has been used to begin the minimization process 
(see Fig. El). 

The Fig. [5] presents the spectral weight function at the unitary regime, for the momentum close to the Fermi 
surface, at the temperature which is slightly smaller than the critical value: T/sf = 0.12. This example represents 
the most difficult situation (from the numerical point of view), where the weight function exhibits the well defined 
structure of two maxima around the Fermi level. However, as one can see, the calculations are surprisingly stable and 
starting from two completely different models one gets practically the same results. It has to be emphasized that the 
results exhibit also the stability of the widths of the spectral weight function. 




FIG. 5: (Color online) Upper subfigures: The extracted spectral weight function (line) at the unitary regime for the momentum 
at the vicinity of the Fermi surface (occupation probability n(p) w 0.6) at T/sf = 0.12 (T c /ef = 0.15(1)). Two sets of results 
have been obtained using MEM for two different models (dashed line) being one gaussian (left figure) and two gaussians (right 
figure). The lower subfigures show the Green's function (open circles) obtained from PIMC calculations and Green's function 
reconstructed from the spectral weight function (line) presented in the upper subfigures. 

The second approach is based on the singular value decomposition of the integral kernel (j3]) [32]. It was shortly 
described in our previous paper [ll[ . Here we mention that the numerical realization of the S VD approach requires to 
introduce two parameters (a and b) which define the interval (a < uj < 6), where A(p,u) is assumed to be localized: 

7 exp(-cjr) f b exp(-cjr) 

duoA(p,uo)- ^— — '-—= / duA(p,uj)- ^— — '—. (7) 

VjP ' j 1 + exp(-cj/3) J a VjP ' j 1 + exp(-cj/3) v J 



f 

J — c 



The prior information concerning the support of the spectral weight function is essential to produce the results of 
high quality [33j| . Since the optimal values of parameters a and b are unknown we have applied the so called bootstrap 
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strategy [34[ . The procedure consists of many reconstructions of the spectral weight function performed for randomly 
generated parameters a and b (within some reasonably chosen interval) . These partial solutions are subsequently used 
to evaluate the final solution (as an average) as well as uncertainties with respect to algorithm parameters. 

Since MEM and SVD methods are based on completely different approaches their agreement serves as a robust 
test for the correctness of the determination of the spectral weight function. In the Fig. [6] the quasiparticle energies 
obtained using MEM and SVD methods have been presented and the quantitative agreement between both methods 
is clearly visible. 
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FIG. 6: (Color online) Quasiparticle energies (maxima of the spetral weight function at fixed momentum) obtained using 
MEM and SVD methods for the unitary limit (left panel) and the BEC side (right panel) at the critical temperature. The 
error bars for the SVD solutions are related to the sensitivity of the maxima localization with respect to the parameters of the 
SVD algorithm and Monte Carlo uncertainties. In both cases the presence of the gap is revealed. Dashed line denotes the fit 
of the BCS dispersion relation to the extracted quasiparticle energies. 



We have performed a number of tests to quantify the limitations associated with MEM and SVD approaches. In 
both cases we have identified the resolution limit which impose the lower bounds on the value of the gap which can 
be detected by both methods. In the case of cold atomic gases in the unitary regime the resolution limit is about 

0. 196f which turns out to be sufficient to reveal the existence of the pseudogap phase. It means that both methods 
provide in practice only a lower bound on the temperature T*. 

Similarly like in Ref. [11] we have fitted the formula for the quasiparticle energy E(p) to the maxima of the spectral 
weight function. The extracted values of gap A, self-energy U and effective mass m* have been presented in the table 

1. The comparison of the results provided in the Table I and those obtained in Ref. [ll| for a smaller lattice allows to 
estimate errors related to the finite size effects. The single-particle parameters do not vary more than by 10%, except 
for the gap at the temperatures close to T*, where it is on the verge or below the resolution limit. The resulting 
uncertainty of the value of T* have been included in the error bars presented in the Fig. 3 in the article. 



ENERGY DISTRIBUTION CURVES FROM THE SPECTRAL WEIGHT FUNCTION 



Since the experiments with cold atoms are performed in a trap one has to translate the results obtained for the 
uniform system to the nonuniform one determined by the geometry of a trap. The quantity of interest is the energy 
distribution curve (EDC), which is measured by exciting atoms from strongly interacting to weakly interacting states. 
Hence the results depends on the occupation of the interacting state and thus are related to the spectral weight 
function: 

EDCfe 25, T) = Cp 2 l°° dr r 2 -^A [^-, E - p.{r) T_\ f{E _ (r)) (g) 
Jo £f(t) Lp F (r) e F (r) e F (r)l 

where /cf(^), /i(r), £f{t) are the local Fermi momentum, chemical potential and Fermi energy in the trap, respectively. 
Here the spectral weight function is expressed in units [e^ 1 ] and parameter C is chosen to make the total area of the 
EDC signal equal to unity, and f(E) = 1/(1 + e@ E ) is the Fermi distribution function at the temperature T = 1//3 
which acts as a filter for the occupied branch of the spectral function. 
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In order to calculate the above expression one has to determine the profile of atomic density in a trap. This has 
been done using the local density approximation (LDA) [25]. In this approach, the grand canonical thermodynamic 
potential for a system confined by a trap potential U(r) is a functional of the local density n(r) given by 



dV 



-Sf(y)(p{x-, fc.p(r)a)n(r) + U (r)n(r) — An(r) 



where 



e F (rY 



ft 2 F 3 

£F(r) = 2^ [37r2n(r)]2/3 ' N = 5 £f ^ x > M>M> 



(9) 



(10) 



where F/N is the free energy per particle. The overall chemical potential A and the temperature T are constant 
throughout the system. The density profile will depend on the shape of the trap as dictated by SQ/Sn(r) = 0, which 
results in: 



5Q 
Sn(r) 



S(F - XN) 
Sn(r) 



/i(x(r), fcp(r)a) + U(r) — A. 



(ii) 



At a given T and A, equations (fTU|) and (fTTj) completely determine the density profile n(r) in a given trap for a given 
total particle number. 

For a given profile of the atomic cloud the integral (|5J) can be determined. The results presented here concerns the 
EDC in the vicinity of the critical temperature T c . As an example we show 2D plots of the EDC at unitarity and at 
the BEC side corresponding to (kpa) -1 = 0.2 at the center of the trap (see Fig. [7]). Results were generated for the 
trap used in the experiment reported in Ref. [23|. If the system is normal the maxima of the EDC should reproduce 
the E oc p 2 dependence of the dispersion relation. At the lower temperatures, where part of the system is superfluid 
or in the pseudogap phase the behavior of EDC is more complicated as it contains the mixture of various dispersion 
relations. One can relate the deviation from the p 2 behavior to the presence of other phases in the system. The 
behavior of EDC at high p corresponds to the occupied branch for momenta above the Fermi level. At unitarity 
this behavior can be attributed to the relation n(p) oc 1/p 4 at large p 15, 35|. It is worth noting that the recent radio 
frequency spectra analysis of a trapped unitary 6 Li also revealed the existence of three distinct phases in a trap below 
T c , similarly like in the upper panel of Fig. [7] [36]. 
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FIG. 7: (Color online) Energy distribution curves obtained from PIMC calculations in the LDA approximation for three 
temperatures (defined in the center of the trap) for JILA trap 23]. The left column corresponds to the unitary regime, and the 
right one - to the BEC side corresponding to (/c^a) -1 = 0.2 defined at the center of the trap. Ef — fiwofiN) 1 / 3 is the Fermi 
energy of the gas in the harmonic trap defined by the frequency ujq — {uJxUJyUJz) 1 ^ . In the middle subfigures the temperature 
corresponds to the critical temperature. The insets show the density distribution sections. Different colors correspond to 
different phases: superfluid phase - blue, pseudogap phase - green, normal phase - red. 



